{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Propagating the Uncertainties on the Effective Area (Aeff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"border: 1px solid #000; padding: 5px; background: #ddd\">\n",
    "    <b>Warning:</b>\n",
    "\n",
    "When evaluating the systematic uncertainties due to the effective area, the flux(source_name, emin, emax) method available in pyLikelihood must not be used.\n",
    "\n",
    "The systematic uncertainty must be calculated using the spectral parameter values in the output xml files.\n",
    "\n",
    "If N is the normalization of the source, one obtains:\n",
    "\n",
    "ϵ<sub>syst,max(min)</sub> = |N<sub>max(min)</sub> - N<sub>nom</sub>|/N<sub>nom</sub>\n",
    "\n",
    "where N<sub>nom</sub> is the best-fit value obtained without any scaling file while N<sub>max(min)</sub> is the one obtained using the max(min) scaling files.\n",
    "\n",
    "The value of ϵ<sub>syst,max(min)</sub> can be used to calculate the systematic uncertainty on the flux F:\n",
    "\n",
    "F<sub>syst,max(min)</sub> = (1 ± ϵ)<sub>syst,max(min)</sub>F<sub>nom</sub>\n",
    "\n",
    "This is strictly valid if the normalization of the source is fitted only, otherwise parameter correlation should be taken into account.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As for any data analysis, it is important to consider the systematic errors when performing a LAT analysis.\n",
    "\n",
    "An overview of the sources of systematic errors in the LAT can be found on the [caveats page](http://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html). Among the sources of systematic errors related to our imperfect knowledge of the performance of the instrument, the most important one is due to innaccuracies in the effective area.\n",
    "\n",
    "The systematic uncertainty on the effective area is modeled by ε(E), the maximal relative difference to the nominal effective area Aeff<sub>nom</sub>.\n",
    "\n",
    "As explained in the [caveats page](http://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html), the LAT team has estimated ε(E) by performing several consistency checks between data and IRF predictions. ε(E) depends on how the analysis is performed, as shown in the following figure and table:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='http://fermi.gsfc.nasa.gov/ssc/data/analysis/aeffsystforfssc.gif'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table class=\"styled-table\" summary=\"Sample Aeff Scaling Files and Functions\" width=\"100%\">\n",
    "  <tbody>\n",
    "    <tr>\n",
    "      <th class=\"center\" colspan=\"5\" width=\"100%\">Aeff Uncertainty Functions</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "      <th>Selection (see <a href=\"http://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html\">caveats</a> for more information)</th>\n",
    "      <th>Plot</th>\n",
    "      <th>E = 30-100 MeV</th>\n",
    "      <th>E = 100 MeV - 100 GeV</th>\n",
    "      <th>E &gt; 100 MeV</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "      <td>FRONT, BACK, FRONT+BACK, Joint Analyses (w/ edisp)</td>\n",
    "      <td>Red Curve</td>\n",
    "      <td>3% + 14% x (2.0-log(E/MeV))</td>\n",
    "      <td>3%</td>\n",
    "      <td>3% + 12% x (log(E/MeV)-5)</td>\n",
    "    </tr>  \n",
    "    <tr>\n",
    "      <td>PSF and EDISP Types (Individual)</td>\n",
    "      <td>Blue Curve</td>\n",
    "      <td>10% + 20% x (2.0-log(E/MeV))</td>\n",
    "      <td>10%</td>\n",
    "      <td>10% + 10% x (log(E/MeV)-5)</td>\n",
    "    </tr><tr>\n",
    "      <td>FRONT, BACK, FRONT+BACK, Joint Analyses (w/o edisp)</td>\n",
    "      <td>Black Curve</td>\n",
    "      <td>5% + 20% x (2.0-log(E/MeV))</td>\n",
    "      <td>5%</td>\n",
    "      <td>5% + 10% x (log(E/MeV)-5)</td>\n",
    "    </tr>\n",
    "      \n",
    "  </tbody>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One has to keep in mind that ε(E) defines an envelope in which any function is a valid effective area displacement (though we do not expect extremely abrupt changes with energy; going from min to max should not occur within less than 0.5 in log(E)).\n",
    "\n",
    "As a result, for any function f(E) with |f(E)| < ε(E), (1+f(E)) x Aeff<sub>nom</sub> is a possible valid effective area. The choice of f(E) depends on the analysis, as will be explained below (see [bracketing Aeff](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/Aeff_Systematics.html#bracketing))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to propagate the systematic uncertainties on the effective area in a usual LAT source spectral analysis, one has to compare the results obtained with Aeff<sub>nom</sub> with the results obtained with a possible valid effective area Aeff<sub>sys</sub> = (1+f(E)) x Aeff<sub>nom</sub>.\n",
    "\n",
    "While it seems straightforward, one has to keep in mind an important technical subtlety: some sources in the source definition XML file should not be convolved with Aeff<sub>sys</sub>. This is due to the fact that a usual source spectral analysis uses some information that has been derived using Aeff<sub>nom</sub>. This includes:\n",
    "\n",
    "* the Galactic diffuse model\n",
    "* the isotropic diffuse model\n",
    "* source spectral parameters in a LAT source catalog"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, to perform a fully consistent spectral analysis with Aeff<sub>sys</sub>, one would have to rederive all of the information about the above sources (i.e. derive new Galactic and isotropic diffuse models and a new source catalog). This is computationally prohibitive and, fortunately, there is a simple way to overcome this difficulty.\n",
    "\n",
    "Let's consider, for instance, S, a catalog source in the source model, which spectral parameters are held fixed in the fit to the catalog values. These spectral parameters were derived in the catalog analysis; that is to say that they were the result of a spectral analysis of a ROI containing this source S.\n",
    "\n",
    "Because the catalog analysis used Aeff<sub>nom</sub>, the spectral parameters are such that the convolution of the source spectrum with Aeff<sub>nom</sub> predicts the correct number of counts due to source S. As a consequence, using these spectral parameters with a different Aeff would predict a wrong number of counts.\n",
    "\n",
    "Since the goal of fixing the spectral parameters of S to the catalog ones is to ensure the prediction of the correct number of counts from S, one can see that one should always use Aeff<sub>nom</sub> for this source S, regardless if one is assessing the effect of Aeff<sub>sys</sub>.\n",
    "\n",
    "This reasoning holds for the case in which only the normalization parameter of source S is free in the fit. Because f(E) can depend on E, it induces a change in the spectral shape that can not be absorbed by the freedom of the normalization parameter in order to predict the correct number of counts as a function of energy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Practically, when performing a spectral analysis of a ROI with a non nominal effective area Aeff<sub>sys</sub>, one has to use Aeff<sub>sys</sub> with the sources for which all spectral parameters are free and Aeff<sub>nom</sub> for all the other sources in the source model (including the Galactic and isotropic diffuse models because they contain spectral shape information that can not be modified).\n",
    "\n",
    "Starting with the [Fermitools](http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/), a mechanism for evaluating the Aeff systematics has been implemented which allows the user to modulate the effective area for an individual source. It is done by giving the possibility to specify the scaling (1+f(E)) for each source.\n",
    "\n",
    "The scaling is given in a scaling file which is just a text file with two columns: energy in MeV and the relative scaling of the effective area. Because points are interpolated on a log-log grid (i.e., a power-law is used to interpolate between points), it is strongly recommended to use a rather fine log(E) binning (a good rule of thumb is at least 15 bins per decade)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The path of the scaling file is given in the `scaling_file` attribute of the spectrum element in the XML definition of a source, as shown in the following example:\n",
    "\n",
    "```xml\n",
    "<source name=\"Mrk 421\" type=\"PointSource\">\n",
    " <spectrum scaling_file=\"scaling_function.txt\" type=\"PowerLaw2\">\n",
    "  <parameter free=\"1\" max=\"1e+10\" min=\"0\" name=\"Integral\" scale=\"1e-07\" value=\"1.3\" />\n",
    "  <parameter free=\"1\" max=\"-1\" min=\"-5\" name=\"Index\" scale=\"1\" value=\"-1.6\" /> \n",
    "  <parameter free=\"0\" max=\"200000\" min=\"20\" name=\"LowerLimit\" scale=\"1\" value=\"100\" />\n",
    "  <parameter free=\"0\" max=\"200000\" min=\"20\" name=\"UpperLimit\" scale=\"1\" value=\"100000\" />\n",
    " </spectrum>\n",
    " <spatialModel type=\"SkyDirFunction\">\n",
    "  <parameter free=\"0\" max=\"360\" min=\"-360\" name=\"RA\" scale=\"1\" value=\"166.1138\" />\n",
    "  <parameter free=\"0\" max=\"90\" min=\"-90\" name=\"DEC\" scale=\"1\" value=\"38.2088\" />\n",
    " </spatialModel>\n",
    "</source>\n",
    "```\n",
    "\n",
    "If no scaling file is defined in the XML definition, the nominal Aeff is used by default."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bracketing Aeff\n",
    "\n",
    "Because any reasonably smooth function f(E) with |f(E)| < ε(E) can be used to define a possible valid effective area Aeff<sub>sys</sub> = (1+f(E)) x Aeff<sub>nom</sub>, the user has to choose a set of functions f(E) to assess the systematic uncertainties.\n",
    "\n",
    "The LAT team recommends to use the so-called *bracketing Aeff* method. It consists, for any given observable (flux, spectral index or cutoff, etc...), the 2 functions f(E) that maximize the negative and positive variations of this quantity.\n",
    "\n",
    "The exact choice of the functions is left to the user since it strongly depends on the observables that are measured.\n",
    "\n",
    "For instance, in order to estimate the systematic uncertainty on a spectral break or cutoff, one should choose functions that flip from -ε(E) to ε(E) or from ε(E) to -ε(E) around the measured value of the cutoff or break."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most extreme cases in terms of total Aeff variations (but not in terms of spectral variations) are the ones in which f(E) correspond to + and - ε(E). As examples we provide the scaling files for the extreme cases:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table class=\"styled-table\" summary=\"Sample Aeff Scaling Files and Functions\" width=\"100%\">\n",
    "  <tbody>\n",
    "    <tr>\n",
    "      <th class=\"center\" colspan=\"2\" width=\"100%\">Example Scaling Files</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "      <th>Selection (see <a href=\"http://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html\">caveats</a> for more information)</th>\n",
    "      <th>File Function</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "      <td>FRONT, BACK, FRONT+BACK, Joint Analyses (w/ edisp)</td>\n",
    "      <td><a href=\"p302_aeff_syst_min.txt\">min</a> <a href=\"p302_aeff_syst_max.txt\">max</a></td>\n",
    "    </tr>  \n",
    "    <tr>\n",
    "      <td>FRONT, BACK, FRONT+BACK, Joint Analyses (w/o edisp)</td>\n",
    "      <td><a href=\"p302_aeff_syst_noedisp_min.txt\">min</a> <a href=\"p302_aeff_syst_noedisp_max.txt\">max</a></td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "      <td>PSF and EDISP Types (Individual)</td>\n",
    "            <td><a href=\"p302_aeff_syst_psf_edisp_type_min.txt\">min</a> <a href=\"p302_aeff_syst_psf_edisp_type_max.txt\">max</a></td>\n",
    "    </tr>  \n",
    "  </tbody>\n",
    "</table>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
